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ABSTRACT 

Transit timing variations - deviations from strict periodicity between successive passages of a tran- 
siting planet - can be used to probe the structure and dynamics of multiple-planet systems. In this 
paper, we examine prospects for numerically solving the so-called inverse problem, the determination 
of the orbital elements of a perturbing body from the transit timing variations it induces. We assume 
that the planetary systems under examination have a limited number of Doppler velocity measure- 
ments, and show that a more extensive radial velocity characterization with precision comparable to 
the semiamplitude of the perturber may remove degeneracies in the solution. We examine several 
configurations of interest, including (1) a prototypical non-resonant system, modeled after HD40307 b 
and c, which contains multiple super-Earth mass planets, (2) a hypothetical system containing a tran- 
siting giant planet with a terrestrial-mass companion trapped in low-order mean motion resonance, 
and (3) the HAT-P-13 system, in which forced precession by an outer perturbing body that is well 
characterized by Doppler radial velocity measurements can give insight into the interior structure of 
a perturbing planet, and for which the determination of mutual inclination between the transiting 
planet and its perturber is a key issue. 
Subject headings: Extrasolar Planets, Data Analysis and Techniques 



1. INTRODUCTION 



While the overall census of extrasolar planets con- 
tinues to climb steadily (453 as of this writingM), the 
emerging population of Earth and Super-Eartn sized 
{Ai sini < lOA^©) planetary companions that has been 
uncovered by high-precision radial velocity (RV) sur 



variations (TTV) will be caused by gravitational pertur 
bations exerted by additional planets, ca using deviations 
from strictly p eriodic Keplerian orbits ( Miralda-Escude 



2002 Holman & Murray 2005). These can be used to 



inter the or bital elements ot the perturbing planet (Agol 



veys (e.g. fRivera et al.||2005l |Udry et al.]|2007[ [Mayor 



i 



et al. 2005), or at least place limits on the presence of 



et al. 2009a " Vogt et al. 20fop is shifting the interest ot 



additional planets (e.g. Alonso et al. 2008 Miller-Ricci 



et al. 2008 1 . An approximate analytic estimate tor I'l'V 



many planet search programs towards terrestrial planets. 
Future refinements in ground-based RV programs will 
likely continue to further push the detection capabilities 
towards the low-mass end of the planet ary population 



amplitude for a tr ansiting planet and an e xternal per- 
turber is given by (Holman fc Murray||2005 ) 



(Mayor et al. 



2009b Howard et al. 



2010) 
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On the other hand, the availability of ground and 
space-based surveys dedicated to photometric monitor- 
ing of large samples of host stars is affording constraints 
on the true mass and bulk composition of the Super- 
Earth planetary population ( Leger ct al. 2009 [Queloz 

et al.|[2009 Charbonneau et aT] 2009). in pa rticular, the 
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Kepler mission (e.g. Koch et ai.||2004 
to yield transiting Earth-mass planets m t 
Zone (HZ) as part of its mission objectives, through con- 
tinuous and simultaneous photometric sampling of more 
than 100,000 dwarf stars. However, this class of objects 
will likely represent a small percentage of the detections 
(given the constraints of the mission design) , and a large 
number of Neptun e-mass and giant plane ts will be de- 
tected as well (e.g. Borucki et al.||2010a|b ). 

The exquisite precision and sheer size of the Kepler 
transit timing datasets of giant planets, as observed dur- 
ing the projected four years to six year mission dura- 
tion, opens up an alternative route to the detection of 
low-mass planetary companions. Indeed, transit timing 



(where we use the symbols A4 for mass, P for pe- 
riod, e for eccentricity and a for semi-major axis; ae = 

^trans/ [(^pert\^ ^pert)\)- 

The amplitude of these variations can be quite large 
and amenable to detection, eit her in the presence of 
high-eccentricity perturbers (e.g. Steffen & Agol 2005 ) or 
when the two planets lie near a low-order mean motion 
resonance (MMR). Indeed, MMRs are an entirely plau- 
sible outcome of core-accretion models of planetary for- 
mation, whereby planets can be captured and locked into 
an MMR during the mig ration stage (e.g. Nelson fc Pa- 



paloizoul 2002 IPapaloizou fc Szuszkiewicz„2005[ [Beauge 
et al. 2006p . Observationally, several of the detected ex- 
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^ exoplanet.eu, retrieved on May 12, 2010 



trasolar systems with multiple planets may be locked 
in low-order MMRs. Three such systems (HD82943, 
HD73526 and HD128311) are engaging in deep 2:1 reso- 
nances well characterized by the observations, and GJ876 
has recently been reported as a Laplace-type resonance 
chain (4:2:1; [Rivera et aH2010| ). For instance, the TTV 
amplitude induced by an Earth-mass perturber in a 2:1 
resonance with a 3-day Jupiter-mass pl anet, both in cir- 
cular orbits, is of order of minutes ( jAgol et al. 2005) . 
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This is a large signal compared to an accuracy i n the 
measureraent o f the central transit time of order (Ford 
fc Gaudi|[2006t 
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(where ig is the duration of the transit ingress/egress, F 
is the observation rate, aph is the photometric precision, 
Rpi and i?* are the radius of the planet and the radius 
of the star, respectively), amounting to 10s of seconds 
for milli-mag photometric accuracy. The recently pub- 
lished JtTepfer central times ([Latham et al .'2009'; Borucki 
et al ||2 010b[ [ Jenkins et al.f|i>0i01 IHuHEam ct al. 2010- 
Koch et al. 2010|) are in rough accordance with this esti- 
mate. _b'urtlicrmore, with respect to the Kepler project, 
we note that once a transit is detected with sufficient 
signal-to-noise ratio, the star will be switched from the 
long-cadence ( 30 minute) to short -cadence (1 minute) 
sampling rate ( Borucki et al.|[2008 ), improving the tem- 
poral resolution ot the transit even further. We take 
<^tr,K = 2 X 10^'' d (sa 15 s) as a conservative estimate 
of accuracy on the central transits. 

Given a large dataset comprising 1 year or more of con- 
tinuous transit monitoring, is it possible to infer the mass 
and elements of the perturbing planet? Reconstructing 
the properties of the perturber from a noisy T TV signal 



is a complex, and p ossibly highly degenerate (Nesvorny 



& Morbidelli 2008), inverse problem. In this paper, we 
present a series of simulations aimed at detecting low- 
mass perturbers from realistic central transit and follow- 
up RV data. To this end, we produce a large sample of 
Kepler-like observations and attempt to characterize the 
perturber using the algorithm tool set offered by a revised 



version of the Systemic Console ( Meschiari et al. 2009 



hereafter Paper I). A number of ditterent planetary re 
alizations were used, in an attempt to fully capture the 
complexity of TTV fitting, drawing the orbital elements 
from observed planetary systems. For the sake of sim- 
plicity, we focus on two-planet systems, but the method 
is fully general within the constraints of CPU time and 
measurement errors. 

The plan of the paper is as follows. In ^ we briefly re- 
view describe the algorithms used to derive best-fit mod- 
els and accompanying error estimates. In f|3] we examine 
the characterization of planets similar to HD40307c and 
d, which lie close but not quite in a 2: 1 MMR. Our analy- 
sis makes use of the HARPS dataset ( Mayor et al.|2009a l 
and a simulated transit timing dataset. In ^we ht the 
synthetic realization of a planetary system deep in a 2:1 
MMR, with an external perturber (using IIAT-P-7 as 
our model system). Finally, in qS^we analyze constraints 
placed by TTVs on the three-dimensional configuration 
of planetary systems, using HAT-P-13 as a test case, and 
conclude in f|6l 



2. NUMERICAL SETUP 

The transit timing variation signal is defined as the 
difference between the observed central transit times and 
the predicted times from a linear regression (correspond- 
ing to a single-planet Keplerian fit with period Pi): 



The variations originate by the mutual gravitational in- 
teractions with an additional body, chiefly causing short- 
term oscillations wherein the true anomaly /i trails or 
leads the Keplerian v alue and long-term effe cts such as 
pericenter precession ( Heyl fc Gladman|2007 ) . In princi- 
ple, since the signal will depend on the INewtonian evo- 
lution of the planetary system, TTVs can provide a sen- 
sitive probe for the three-dimensional orbit of the second 
planet, in combination with the tight constraints on the 
eclipsing planet's period and the time of pericenter pas- 
sage provided by the central transits themselves. 

However, solving the inverse problem of deriving a 
best-model fit to the TTV observations can be daunting. 
The computation of the predicted TTV signal requires 
precise N-body integrations (with N > 3). In the gen- 
eral case, the dependence of the signal on the set of or- 
bital parameters is not directly clear; unlike, e.g. the RV 
technique, deviations from the Keplerian signal - as op- 
posed to the Keplerian signal itself - constitute the bulk 
of the information. The use of Fourier analysis to sort 
out periodicities in the data is generally hampered by 
the sparseness of the transit observations. Furthermore, 
given the extreme sensitivity of St to the model parame- 
ters, local minimization routines can easily get stuck in 
narrow x^ minima, or fail due to steep gradients in the 
landscape. Finally, as shown in the later sections, there 
is a degree of non-uniqueness as multiple models can fit 
the transit timing observations when mea surement er- 
rors are ta ken into account (see also e.g. Nesvorny & 
Morbidelli 2008) ; these degenerate solutions are charac- 
terized by comparable x^ ^ 1? ^-nd must be taken into 
account when deriving parameter uncertainties. 



Direct se arches of the parameter space (e.g. Steffen & 
Agol||20(y7 ) can be extremely expensive in ternis ot CP U 
time. A more appealing alternative is represented by the 
TTV Inversion Method (TTVIM; [Nesvorny fc Beauge 
2010 and related papers). TTVIM combines a fast algo- 
rithm for computing the 2-planet transit timing based on 
perturbation methods with a downhill simplex method to 
obtain good convergence towards the perturbing planet's 
parameters. However, some issues remain in addressing 
systems lying close to a MMR. 

In this paper, we adopt the approach of finding best- 
fit models to joint TTV and Doppler velocity data sets 
by driving an efficient Bulirsch-Stoer integrator with the 
Simulated Annealing algorithm integrated in the Sys- 
temic Console (Paper Irj SA-type algorithms are well- 
suited to exploring the orbital parameter space (period 
P, mass M, eccentricity e, inclination i, mean anomaly 
at epoch Mq , longitude of pericenter zu and node f2 for 
each planet) and converging, in principle, to global min- 
ima (subject to appropriate choices of scheduling algo- 
rithm and scale parameters). Several minimizers can be 
run in parallel with different initial temperatures and ini- 
tial conditions, exploiting modern multi-core CPUs capa- 
bilities. The step size vector is automatically adjusted to 
attain an acceptance rate of ~ 25%; we have found em- 
pirically that this value is an optimal compromise. After 
a fixed number of ste ps, we invoke a do wnhill simplex 
algorithm (AMOEBA; [Press et ar][1992|) in an attempt 



Stk 



tk - kPi 
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^ The new version of the Systemic Console, including a Bulirsch- 
Stoer integrator, AMOEBA and fully non-coplanar fitting is avail- 
able for download at www.oklo.org 
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Fig. 1. — Sensitivity of the RV method to the mutual gravi- 
tational perturbations: Keplerian model subtracted from the in- 
tegrated model (thick curve) compared to the HARPS residuals 
(empty circles). 

to home in on nearby deep minima. This avoids missing 
promising solutions when the SA step size is too large 
to properly resolve them. In practice, this scheme per- 
mits the derivation of the full set of degenerate solutions 
compatible with the observational errors. 

Although we recognize that this approach can be com- 
putationally inefficient compared to perturbation meth- 
ods, the implementation is trivial and can use existing 
integration techniques. Furthermore, it permits the char- 
acterization of arbitrary planetary configurations (in- 
cluding Npi > 2, resonant, high-eccentricity and inclined 
bodies) and the inclusion of additional dynamics (such 
as tidal evolution) self-consistently, owing to the fully 
general N-body integration. Finally, we remark that in 
this work the parameters of the transiting planet are not 
fixed, but derived simultaneously from the available data. 
This mimics follow-ups of transiting planets, whereby the 
mass of the transiting planet is determined by a small 
number of RV measurements. 

We use the combined x^ statistic detailed in Paper 
I to simultaneously fit the transit timing and follow-up 
RV datasets. While there is a degree of ambiguity in the 
choice of the weighing factor A, this is not a concern in 
the vicinity of a solution, where the contribution from 
RVs and transits is approximately equal for A = 1. Far 
from the solution, the contribution from transits to the 
X^ budget is extremely large; however, this is not an issue 
in practice because we first fit for a one-planet solution, 
reducing the initial x^ to x^ ~ 5t/(jTR- 

3. HD40307 



Mayor et al. 



(2009a) recently announced a three 
Super-Earth planetary system orbiting the nearby metal- 
deficient dwarf HD40307. Interestingly, while this system 
lies close to a 4:2:1 Laplace resonance chain, such a con- 
figuration is ruled out by the observations. The a-priori 
transit probability for the innermost 4.3d planet is high 
enough to warrant a transit follow-u p; unfortunately, n o 
transit was detected using Spitzer (Gillon et al. 2010), 
preventing the placement of desirable constraints on the 
bulk composition of the three planets. 

We derived the orbital parameters for the system using 
the publicly available HARPS dataset, obtaining best- 
fit and error estimations in good accordance with the 
published configuration. As can be seen by Figure [T 
the difference in RV signal between a fully integratec 
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Fig. 2. — Predicted transit timing variations for planets b (empty 
circles) and c (black circles), over the HARPS observation window. 

model (using Bulirsch-Stoer) versus simple superposition 
of Keplerian orbits is negligible compared to the HARPS 
error bars and the RV residuals. 

As a comparison, we derived synthetic TTV computed 
comparing a simple linear fit to synthetic transits com- 
puted with the fully integrated solution above; we as- 
sumed, respectively, planets 6 and c to be transiting and 
computed the primary transit timings for the HARPS ob- 
servation window. To each transit timing observation, we 
added a Gaussian white noise of amplitude a = 2 x 10"'* 
d (0.3 minutes), as a simple, conservative model for Ke- 
pler timing uncertainties. The TTV dataset is shown in 
Figure [2] The amplitude of the TTV signal for planet 
c is approximately Scr, making it a far more sensitive 
probe of the mutual gravitational perturbations than the 
highest-precision RV measurements available. 

Although no transits have been so far detected for 
planet b, the transit probability for planet c is a tantaliz- 
ing 5% and the transit depth is of order 400 ppm, fully 
within the capabilities of Kepler. Therefore, it is an in- 
teresting illustrative test-case problem to use the known 
orbital elements of the HD40307 system and analyze the 
constraints imposed by TTV on the perturbing planet, 
in absence of high-precision radial velocities. Given that 
the bulk of the signal originates from the mutual per- 
turbation between planets c and d, we hereafter solve 
the simpler two-planet inverse problem and neglect the 
contribution from planet b. The orbital elements of the 
generating fit are reported in Table IT] 

We generated two sets of central transit observations 
spanning 100 days (11 transits) and 365 (38 transits); we 
assumed every transit is detected with atr = 2 x 10~^ d. 
We also computed a small set of "follow-up" synthetic ra- 
dial velocity observations (10 points), which set the scale 
for the mass and the eccentricity of the transiting planet 
(period and mean anomaly at epoch being primarily de- 
termined by the transit timing). We draw the measure- 
ment errors to mimic mid-range precision observations; 
the average measurement error is ^ 1.5 m/s. As a com- 
parison, we also computed a third synthetic RV dataset 
drawing from the HARPS schedule and measurement er- 
rors for the two planets alone; we assumed a small jitter 
of ~ 0.7 m/s. We point out that this jitter is excellent, 
and depending on the properties of the parent star, a 
realistic case might require more RVs. 

Eight SA simulations were launched (one per core on 
a Mac Pro Xeon workstation), with initial temperature 
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TABLE 1 
Best-fit solutions 



Best fit (HARPS) Best fit (lOOd) Best fit (365d) 



P (days) 
M {Mj) 
e 



9.621 (1) 

20.439 (5) 

0.0218 (6) 

0.0290 (8) 

0.06 (3) 

0.12 (2) 

284 (6) 

12(7) 



9.6214 (4) 
20.2 (2) 
0.021 (5) 
0.025 (4) 
0.036 (3) 
0.06 (3) 
358 (4) 
78 (23) 



9.62114 (5) 
20.45 (1) 
0.02 (1) 
0.025 (2) 
0.034 (4) 
0.01 (2) 
358.2 (1) 
71(4) 



X^ 






10.49 


1.29 


1.15 


RMS 


(m s 


-') 


1.04 


1.17 


1.25 


Xtr 






— 


0.4 


0.75 



Mark ov Chain Monte Carlo algorithm (MCMC; |Ford 
2005| . While the synthetic HARPS dataset is amenable 
to a bootstrap resampling technique, the rugged x^ land- 
scape for the TTV dataset turned out to be excessively 
complicated for an efficient exploration, yielding artifi- 
cially low parameter uncertainties. The simple MCMC 
algorithm presented in Paper I derived error bars in ac- 
cordance with bootstrap estimates for the RV dataset. 
We use uniform priors in {logP, logM, Mq, e, lu}; while 
more sophisticated approaches are available (e.g. incor- 
porating information from Eqn. [lias a constraints), the 
size and precision of the synthetic aatasets provide strong 
constraints on the model parameters and the choice of 
the priors should not affect our results (Ford 2006). We 
construct MCMC chains of 50,000 states, each state con- 
sisting of 200 iterations. The initial 10% portion of the 
chain is considered "burn-in" and discarded. 
The best-fit solutions to the three datasets (synthetic 
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Note. — Best-fit solutions for tlie HD40307 system. Tfie error on 
the least significant digit is indicated in parentheses. 

and step size regulated to achieve 25% acceptance rate 
in each orbital parameter. The initial configuration used 
the parameters from a single planet best-fit for the tran- 
siting planet, and random elements for the perturbing 
planet (period, mass, eccentricity, mean anomaly and 
longitude of pericenter), avoiding orbit-crossing config- 
urations. For the sake of efficiency, we constrained the 
period of the second planet between 1.5 and 5 times the 
period of the inner planet, the masses between 0.3 and 32 
A^0 and the eccentricity between circular and 0.5. This 
parameter range approximately spans the region where 
the transit timings are sensitive to the perturbations, 
but the reflex stellar semiamplitude K is not so large 
to be readily picked up by RV observations. Finally, ev- 
ery 2,000 steps the current configuration was submitted 
to the AMOEBA routine to attempt direct convergence 
to a solution. The minimization routine is considered to 
be converged and the solution is retained if Xt_r < 1-1 
and the radial velocity RMS < orv, corresponding to 
a combined x^ ~ 1.3. After a predetermined number 
of steps (10,000), if no improvement in the total x^ has 
been reached, the suboptimal solution is discarded and 
a new set of initial conditions is chosen. A sample of 20 
candidate solutions was derived for each dataset; each 
solution representing a local minimum. The lowest x^ 
solution was chosen as the representative best-fit. 

An estimation of the uncertainties on the orbital pa- 
rameters of the best-fit solutions was derived usin g the 
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Fig. 3. — Results of a MCMC simulation consisting of 50,000 
states computed from the HARPS dataset (green points), 100 days 
of transit timing observations -I- 10 follow-up RVs (red points), 
365 days of transit observations -|- 10 follow-up RVs (blue points). 
The parameters of the originating system are marked with a star 
symbol. 



HARPS, 100-d and 365d TTVs + RV foUowup) and re- 
spective uncertainties are compared in Table [l] To our 
knowledge, this is the first attempt to fully fit and derive 
error estimates on a large TTV dataset. We show the 
parameter scatter for the second planet in Figure [3] 

The computed parameter uncertainties show a number 
of interesting properties. Firstly, the period and mass of 
the second planet are derived to an accuracy compara- 
ble to that of the full HARPS dataset, which spans 4.5 
years. The detection of a low-mass planet at this level of 
accuracy showcases the potential of scanning the future 
Kepler datasets for TTV detection candidates. Once 
again, we stress that our estimate of the central transit 
timing noise is likely conservative and that stars on the 
short-cadence list will be observed with an even higher 
accuracy. While the period of the transiting planet is 
constrained by the transits timing themselves, the mass 
is not well constrained because, to a good approxima- 
tion, the amplitude of the TTVs does not depend on the 
mass of the transiting planet itself (Equation [I]) in the 
non-resonant regime. 

Finally, we remark that although the x^ landscape al- 
lowed for several, well-separated local minima, both the 
SA and the MCMC algorithms were able to efficiently 
sample the parameter space. Therefore, it is likely that 
global minimization routines will be part of the standard 
toolset to analyze the future Kepler transit datasets. 
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Fig. 4. — Libration (in degrees) of the resonant arguments ©i 
and Aii7 for the reference configuration (black line) and the best- 
fit configuration (grey line). 

4. HAT-P-7 

The bright nearby dwarf HAT-P-7 hosts a transiting 
hot Jupiter first characterized by the HATNet project 



(Pal et aL 2008). The star is in the field of view of 



one of the Kepler detectors; ten days of photometric 
data, as processed by the Kepler p ipeline, were obtaine d 



during the commissioning phase (Borucki et al. 2009). 



Additional primary transits and a number of sec ondary 



eclipses were observ ed using EPOXI and Spitzer ( Chris 



tiansen et al.||2010[ ), with the intent of studying the at 
mospheric properties of the planet. The EPOXI best-fit 
central times achieved an accuracy oiatr ~ 10^'^ d (w 1.5 
minutes) . 

Given its extensive and diverse coverage, and the in- 
clusion of this planet in the Kepler star list, we chose this 
system as a prototype of the class of massive transiting 
planets that will be monitored by the Kepler mission and 
may reveal TTVs. In particular, we are interested in as- 
sessing the secure detection of a low-mass planet in a 2:1 
MMR with the transiting gas giant (we consider only the 
case of an external perturber in the present analysis) . 

We generated a realistic resonant configuration self- 
consistently with the following procedure. We placed the 
two planets (denoted as 1 and 2, respectively the tran- 
siting planet and the external per turber) on or i ginall y 
widely separated orbits; following Lee & Peale (2002), 

=4^ 



we added a forced migration {a/a = —3 x 10~' yr" 
and an eccentricity damping (e/e = 100 a/ a) term of the 
outer planet to the equations of motion until resonant 
capture is achieved. In this reference configuration, the 
outer planet was captured into an antialigned configura- 
tion with 01 = 2A2 — Al — n72 librating around 0° and 
Atn = ■W2 — TJJi librating around 180°, with an ampli- 
tude of « 5° (Figure [4]). The final eccentricities for this 
choice of forced migration terms are low (ei — 0.002, 62 
= 0.027). 

To illustrate the process, we chose a mass for the sec- 
ond planet of « lOM^ , since this can yield a TTV signal 
larger than 1 minute, easily detectable with Kepler. Fig- 
ure [5] shows the amplitude of the TTV signal for a choice 
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Fig. 5. — (Top) Grayscale map of St (in units of cFtr,K ~ 15 
s) for 10,000 realizations spanning a range of perturber periods 
and masses, using the reference configuration for the other ele- 
ments. (Below) As above, in the region near the 2:1 resonance. 
The contours show the parameter space where St > 2crtr, assum- 
ing atr = 2 X lO^** d {Kepler, red contour) and atr = 10""^ d 
{EPOXI , blue contour) respectively. The star symbol represents 
the reference configuration. 
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Fig. 6. — Best-fit solutions for the HAT-P-7 dataset lying near 
the 2:1 resonance (circles) and the 3:1 resonance (squares), for two 
different levels of noise in the TTV measurements: 2x 10~* (empty 
symbols) and 5 X 10~^ (filled symbols). 
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of periods and masses, at fixed eccentricies and phases; 
as expected, the TTVs are largest in the proximity of 
resonances. In particular, 3:2, 2:1 and 3:1 MMRs yield a 
sizable TTV signal for our range of perturber masses. 

We created a TTV dataset spanning 1 year (166 ob- 
servations) following the procedure in Section pi using 
the reference configuration as our generating system and 
Gaussian noise at the level of 2 x lO"'* d. We drew 
from the schedule and uncertainties of t he Keck/HIRES 
follow-up observations (Pal et al. 2008) to generate the 
accompanying RV dataset. We note that given the small 
senn-amplitude K2 (~ 2.8 m/s, larger than the typical 
error in the Keck dataset but smaller than the stellar jit- 
ter sa 3.8 m/s) and the few RV points available, the RV 
dataset places only a weak constraint on the parameters 
of the perturbing planet. 

We launched a number of SA chains and allowed the 
parameters of the perturbing planet to fioat freely. We 
found that the best-fitting solutions comprised a set of 
degenerate configurations, shown in Figure[6j The fitting 
routine found two groups of solutions: configurations ly- 
ing near a 2:1 MMR and configurations lying near a 3:1 
MMR can fit the TTV signal equally well. Addition- 
ally, the degeneracy between mass and eccentricity of the 
perturbing planet makes it impossible to place a strong 
constraint on the mass of the second planet. 

This non-unique ness of the inverse proble m was al- 
ready noted in Ncsvorny & Morbidelli ( |2008 [); the mea- 
surement errors filter out some of the 'i"i'V harmonics. 
The authors also pointed out that the non-uniqueness 
threshold (the measurement uncertainty that leads to a 
unique solution) of the number of transits detected; ac- 
cordingly, we verified that a transit dataset covering 2 
years of observations still yielded the two groups of solu- 
tions. Reducing the error on the transit measurement to 
5 X 10~^ d (4 seconds), while not breaking the resonance 
degeneracies, reduced the range of possible masses some- 
what (Figure J6]) . Finally, only a fraction of the solutions 
(about 10%) have librating resonant arguments; the ones 
that do show a much larger amplitude of libration than 
the reference system (61 - 20 - 40°, Atu - 30 - 70°; 
see Figure I4 1. This suggests that the TTV signal alone 
is not enough to constrain the resonant angles. 

Our result is particularly remarkable in that the best- 
fitting solutions cluster around two different MMRs, pre- 
venting a precise characterization of the resonance. Since 
the two solutions yield a different RV semi-amplitude K2 , 
this degeneracy may be broken with RV observations. 
Even a small RV dataset, where uncertainty and jitter 
do not completely wash out the planetary signal, can 
help constrain the parameters of the perturbing planet 
to a reasonable level. Indeed, we verified that a second 
RV dataset comprising 20 measurements with lower jit- 
ter (~ 1 m/s) sufficed to constrain the best-fit solutions 
to the neighborhood of the 2:1 MMR. We conclude that 
while TTVs can be usefully exploited to infer the pres- 
ence of low-mass perturbing planets, a small number of 
RV measurements with a precision comparable to K2 is 
crucial in recognizing the nature of the planetary com- 
panion. This fact makes it much more desirable to find 
configurations orbiting bright parent stars. 

5. HAT-P-13 




300 320 340 

JD - 2,455,000 



Fig. 7. — Sample TTV signals for four different inclinations of 
HAT-P-13 c. Orbits close to perpendicular give rise to large TTV 
signals. 



HAT-P-13 was the first system known to contain a 
transiting planet, &, and an eccent ric outer plane t, c, 
well characterized through RVs (Bakos et al.||2009[). No 
transits of planets c have been detected thus far. A com- 
plete characterization of the three-dimensional configu- 
ration of th e system can establi sh the internal structure 
of planet h ( Batygin et al.|2009 ) and possibly the forma- 
tion and scattering history of the system, with certain 
ranges of inclina tion being favored on theoretical grounds 
( |Mardhng|[2010 |. 

'iransit timing variations can provide the required con- 
straints on the mutual inclination (/) and the nodal line 
marking the intersection of the two orbital planes (fi), 
should transits of c not be detected. The amplitude and 
shape of the TTV signal depe nd significantly on the two 
parameters ([Payne et al. 2010 1, although this dependence 



IS not triviaE 

Figure [7] shows the TTV signal for a number of inclina- 
tions. We centered our dataset around Tp^ri. c since the 
different solutions can be best distinguished by the sharp 
feature in the neighborhood of the pericenter passage of 
c. While the discovery paper predicted a TTV ampli- 
tude of orde r 15-20 seconds, th e updated configuration 
presented in Winn et al. (2010) reduces the expected St 
near the pericenter p assage by a factor ~ 2, to about 



7 seconds for / w 0. Winn et al. (2010) also measured 



a prograde Rossiter-McLaughlin ettect, suggesting that 
both orbits are prograde. 

We produced several transit datasets for mutual incli- 
nations in the range 0° < / < 90° and r2 = 0, assuming 
that all transits between Tpf^^i, c ~ 100 d and Tp^^ic + 100 
d are detected; the other elements were dra wn using the 
published uncertainties (Winn et al. 2010). We added 
white noise to the TTV signal at the 4 x 10"^ d = 3.5 
s level (in order to have St/atr > 2). The RV measure- 
ments were generated drawing from the schedule and un- 
certainties of the Keck/HIRES dataset as reported in the 
discovery paper. 

We used our usual fitting procedure (Bulirsch-Stoer 
as our integration scheme and Simulated Annealing and 
AMOEBA in tandem to pinpoint the solution) , with the 
published fit as our starting configuration. When a so- 
lution was found, we estimated the uncertainty by run- 
ning our MCMC algorithm. We generated 4 x 10^ trial 
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Fig. 8. — Relative inclination distribution for synthetic HAT-P-13 
realizations with 7 = 0,5, 15, 45 and 75° respectively. The median 
inclination and standard deviation are given inside each plot. 

models; of those, the first 10% was discarded and only 
one model every 50 was retained in order to minimize 
correlations between successive elements of the Markov 
chain. Figure [8] shows the marginal distribution of the 
fitted relative inclinations for systems with various de- 
grees of inclination (/ = 0°, 5°, 15°, 45° and 75°). The 
inclination is well constrained for polar and near-polar 
configurations of the outer planet, where the TTV sig- 
nal is sizable; on the other hand, for low inclinations 
there is a large range of allowed configurations. How- 



ever, it is clear that while the inclination distributions 
are broad, they are consistent with the originating con- 
figuration and can discriminate between low-inclination 
and high-inclination configurations. 

6. DISCUSSION 

In this paper we outlined a procedure to solve the in- 
verse problem of deriving best-fitting model parameters 
and associated uncertainties using synthetic radial veloc- 
ity and transit timing variations datasets simultaneously. 
The procedure exploits a number of numerical algorithms 
that are made available to the community through the 
Systemic Console package. 

We tested our fitting method against a number of syn- 
thetic realizations of different planetary configurations, 
including a system of non-resonating coplanar super- 
Earths, a system in a deep 2:1 resonance and a non- 
coplanar system. The transit timing datasets were de- 
rived assuming continuous photometric coverage as pro- 
vided by Kepler, and thus are fully realistic to the extent 
that the transit timing error can be modeled as white 
noise with a constant amplitude. Our analysis shows 
that combined RV and TTV datasets carry enough dy- 
namical information to characterize a system in its full 
three-dimensional configuration. 

Inverse problems have a storied place in astronomy, 
with the discovery of Neptune providing a canonical ex- 
ample. In that case, a fortunate orbital geometry al- 
lowed Neptune's sky position to be pinpointed with suf- 
ficient accuracy that the "prediction" of a new planet 
could credibly be claimed. It is worth pointing out, how- 
ever, that the accurate ephemeris for Neptune in 1846 
was something of a lucky accident. Both Adams' and 
Le Verrier's m asses and semi-major axes were badly off 
(Grant 1852). The correct position of the planet that 
emerged troin the calculations stems from a degeneracy 
of solutions during the period surrounding the conjunc- 
tion of Uranus and Neptune. 

We have found that a similar state of affairs might ap- 
ply to the transit timing measurement scenarios that will 
emerge from Kepler. While departure from strict period- 
icity can be readily measurable, it is generally difficult to 
work out the complete system configuration from transit 
timing measurement alone. We confirmed that the sup- 
pression of TTV harmonics by the transit timing noise 
can lead to severe degeneracies i n the model para meters. 



as first pointed out by Nesvorny fc Morbidelli ( [2008} , 
even when very low levels of timmg error is added to 
the synthetic data. In presence of such degenerate set 
of solutions, however, we have verified that adequate RV 
data can single out the correct orbital configuration. We 
note that other constraints derived by extracting more 
observables from the photometry, such as the t he dura- 
tion of the transits and th eir variations (TDV - Kipping 



et al.|200"9 Kipping|2010 ), may also help remove the de 
generacies in the solution. Including these contributions 
will require more sophisticated modelling approaches. 

Finally, we note that our work did not investigate other 
competing effects that contribute to the TTV signal, 
chiefiy including, but not limited to, light travel time, ex- 
citation of tidal modes in the host star, general relativity 
and the presence of additional planets. Furthermore, the 
investigation of planetary systems with Npi > 2 with the 
methods presented here might be computationally costly 



Meschiari et al. 



due to the large parameter space. 
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